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1.0  INTRODUCTION 


The  minimization  of  tunnel  wall  interferences  has  become  one  of 
the  major  tasks  after  the  introduction  of  ventilated  transonic  tunnels. 

A  variable,  but  uniformly  distributed,  porosity  wall  was  designed  to 
reduce  interferences  at  various  Mach  numbers,  e.  g. ,  the  Aerodynamic 
Wind  Tunnel  (4T)  at  AEDC.  The  recent  requirement  for  an  increase  in 
the  size  of  the  testing  model  to  achieve  higher  Reynolds  number  creates 
severe  interference  which  prohibits  obtaining  useful  data.  In  addition, 
the  axial  gradients  of  interference  may  cause  interference  on  pitching 
moment  for  a  long  model.  By  introducing  an  axially  distributed  poros¬ 
ity  in  the  walls  of  a  slotted  tunnel,  the  elimination  of  pitching  moment 
and  lift  interferences  was  achieved  in  the  experimental  development  of 
walls  for  V/STOL  testing  (Ref.  1).  It  is  necessary  to  search  for  a 
theoretical  optimum  porosity  distribution  for  the  minimization  of  inter¬ 
ference  as  the  guideline  for  an  experimental  program. 

The  first  theoretical  approach  to  the  problem  has  been  carried  out 
in  Ref.  2  to  reduce  the  interference  in  a  two-dimensional  perforated 
tunnel  by  a  gaussian  type  distribution  of  porosity  with  an  approximate 
method.  Specifically,  a  system  of  integral  equations  was  derived  using 
Fourier  transform  and  convolution  theorems  and  then  solved  by  the 
collocation  method  with  a  series  form  representing  the  unknown  func¬ 
tions.  The  selection  of  a  gaussian  distribution  is  strictly  based  on  the 
merits  of  mathematical  simplicity.  The  reduction  of  interference  is 
achieved  (Ref.  2)  by  using  a  simple  singularity  to  represent  the  test 
model.  This  has  been  extended  to  a  finite  chord  airfoil  to  permit  com¬ 
parisons  directly  with  experimental  data  (Ref.  3).  However,  the 
approximate  method  is  limited  to  certain  porosity  distributions.  The 
complete  elimination  of  the  magnitude  and  axial  gradient  of  interference 
requires  a  nongaussian  porosity  distribution.  To  provide  such  a  solu¬ 
tion,  a  numerical  method  for  computing  the  interference  has  been  devel¬ 
oped  to  search  for  an  optimum  configuration  in  the  present  study.  The 
application  of  a  modified  method  of  Block  Cyclic  Reduction  (Ref.  4)  to 
the  lift  interference  computation  is  presented.  The  scheme  to  search 
for  an  optimum  configuration  is  discussed  and  extended  to  a  finite  air¬ 
foil.  The  lift  interference  is  calculated  for  an  NACA  64-series  airfoil 
in  an  optimum  configuration  to  demonstrate  the  achievement  of  mini¬ 
mization  of  interference.  The  effect  of  test  section  length  is  briefly 
examined. 
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2.0  GENERAL  ANALYSIS 


The  lift  interference  in  a  two-dimensional  porous  transonic  tunnel 
is  formulated  for  tunnel  walls  with  varying  porosity  distributions.  The 
optimum  porosity  distribution  may  then  be  obtained  by  judicious  selec¬ 
tion  for  a  given  application. 


2.1  FORMULATION  OF  MATHEMATICAL  PROBLEM 

The  field  equation  of  an  inviscid,  irrotational  fluid  for  subsonic 
flow  in  terms  of  the  perturbation  velocity  potential  <&  in  X-Y  coordinates 
(Fig.  1)  is 


B2  Hi  +  Hi  =  Q 
ax2  3Y2 


(1) 


For  the  boundary  condition  of  the  tunnel,  the  average  mass  flow  is 
assumed  proportional  to  the  pressure  drop  across  the  porous  wall  as 


R(x) 


3$ 

3X 


at  Y  =  ±h 


(2) 


where  R{x)  is  the  empirical  constant,  or  porosity  parameter,  of  the 
porous  wall  and  is  a  function  of  streamwise  location. 


MODEL  IN  A  TUNNEL 
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Figure  1.  Boundary  value  problem  for  tunnel  lift  interference. 
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Within  the  assumptions  of  linerarized  theory,  the  perturbation 
velocity  potential  may  be  divided  into  two  parts  as 

$  =  4>  +  <t>m  (3) 

where  <j>  is  the  interference  potential  caused  by  the  presence  of  tunnel 
walls  and  0m  is  the  disturbance  potential  induced  by  a  model.  The 
linearity  of  the  field  equation  and  boundary  conditions  in  the  normal¬ 
ized  coordinates  x  =  X/j3h,  y  =  Y/h  gives 

111  +  Hi  =  o  (4) 

3x2  3y  2 

and 

R(x)  34>  .  30  _  /r(x)  3<|,m  A  8<t>m\  ..  _  fc-x 

T“  TZ  *  37  "  "  \“3  TiT  ±  3y  /'  Y  “  11 

with  the  upstream  and  downstream  conditions  described  as 

4>(±°°)  =  0  (6) 


The  formulation  is  completed  with  the  set  of  Eqs.  (4),  (5),  and  (6). 
The  finite  difference  method  will  be  used  to  solve  this  system.  An 
efficient  numerical  scheme  is  provided  by  the  modification  of  Block 
Cyclic  Reduction  to  yield  a  solution  of  the  finite  difference  equations. 


2.2  FINITE  DIFFERENCE  EQUATIONS 

To  develop  the  finite  difference  equations,  it  is  assumed  that  the 
interference  potential  <j>  effectively  becomes  zero  at  a  large  finite  dis¬ 
tance  from  the  model  location.  This  distance  will  be  denoted  x*. 

Consider  the  rectangular  region 

_  t-x*  <  X  <  X* 

R  i-1  <  y  <  1 

Let  N  be  any  positive  integer  and  let  k  be  any  nonnegative  integer. 
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Define  M  =  2^.  Let  the  region  R  be  overlaid  with  a  rectangular  net 


with  spacings 


6xi  =  xi+l  "  xi 


1  =  0,  _ _  N-l 

where  the  mesh  points  in  the  x  direction  may  be  distributed  as  desired. 
It  will  be  required  that  xQ  =  -x*  and  x^  =  x*. 


In  the  y  direction,  6y  =  i  and  y^  =  j<5y  j  =  0,  ±1,  ± 

For  notational  convenience,  column  vectors  such  as 


±M. 


u  = 


u. 


will  be  denoted  as 


"N 


u  =  col  (u.^  u2,  .  .  .  .  ,  uN) 


and  any  NxN  tridiagonal  matrix  K  of  the  form 


K  = 


bl  °1 


a2  b2  C2 


aN  bN 


will  be  denoted  by  K  =  (a^  bj,  ci)NxN- 
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Let  the  value  of  the  solution  of  the  finite  difference  equations  at  the 
point  (x^,  yj)  be  denoted  as  <t>[  j  and  let 


^N-l, a) 

1  =  0,  +1, 


(V) 


•  a  •  •  f 


±M. 


3  ^  d)  g  2  ^ 

Using  the  centered  second  difference  approximation  for  — 2-  and  — - 

3x2  ay2 

as  given  in  Ref.  5_for  variable  steps,  the  finite  difference  approxima¬ 
tion  to  Eq.  (4)  on  R  becomes 


(' 


♦i+i.i 


♦i,3 


6xi(6xi+6xi_1) 


sxi-isxi 


;xi-i(5V5xi-iV 


(8) 


d> .  .  ,  .  —  26  .  .  +  <().  .. 

+  i’l+1 _ liJ _  J".1...  =  o 

2Sy  2 


i  =  1, . ,  N-l 

j  =  0,±1,....,  ± (M-l) 


Let  =  2  6yz/^6x^  ^6  x^  +<Sxi_1^J 

bi  =  -2  [i  +  syV  («Wi)] 

-  2iy2/[sx._1  (6x.  *lVl)] 


Since  it  is  required  that  $  (-x*,  y)  =</>  (x*  y)  =  0,  there  results 
<j> n  j  =</>0  j  =  0.  j  =0,  ±1,  ....  ±M.  Then  Eq.  (8)  may  be  written 

♦j+1  +  A*j  +  *j-l  =  0 
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where  A  is  the  matrix 


A  = 


ai)  N-l 


X  N-l 


(9) 


The  procedure  along  the  boundaries  y  =±1  is  as  follows: 


Let  b: 


±_73vx i-*1’  awM 

i  \  3x  i  3y  / 


where  T.  =  3/R(x±)  i  =  1,2,...,  N-l 

On  y  =+l,  the  difference  approximations  given  in  Ref.  5  are  again  used 
to  approximate  Eq.  (5)  resulting  in 


Pi*i+1,M  +  <Ji  *i,M  +  ri*i-l,M 

^i,M+l  ^i.M-1 


+T  . 
1 


26y 


=  B+ 
l 


where  p±  =  jx^/  [ 6x±  ^6x±  + 

q.  =(6x.  -  6xi_1)/(5xi6x._1) 


(10) 


and 


r . 

l 


=  -  5x 


(5xi +  sxi-i)] 


Eq.  (4)  is  also  required  to  hold  for  y  =  +1  which  gives 


i,M+l  +  ai  +1+1.M  +  bi  ♦i.M  +  Ci  ♦i-l.H  +  ♦i.M-l  *  0 


Eliminating  <j>^  from  these  two  equations  yields 

T  *M-1  =  E*M  +  ?+ 


(ID 
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where  E  is  the  matrix 
E 


[2  (r^  -  Tici)  ,  2  {qi  Tibi)  '  2  (pi  Tiai}]  N-l  x  N-l 

with  q*  =  26y  q±,  p*i  =  2<5y  p^  r*  =  26y  ft  =  -5y  Bt 
and  where 

1  -  (°'  Ti'  °)  N-l  X  N-l- 

In  a  similar  manner  it  can  be  shown  that  on  the  boundary  y  =  -1, 


Th-M  *  EJ-M  +  f 


(12) 


The  set  of  finite  difference  equations  (Eqs.  (8),  (11),  and  (12))  is 
readily  solved  for  the  determination  of  lift  interference  once  the  lift 
potential  is  established. 


3.0  LIFT  INTERFERENCE 


The  lift  interference  factor  is  defined  by 

-  C  1 

SCL  U  3y 

In  particular,  the  factor  along  the  centerline,  y  =  0,  can  be  obtained 

by 

,  _  c  1  h  ■  $- 1 

SCL  U  2«y  <13) 


where  <f>  ^  -  6  can  be  computed  by  solving  an  N-l  system  of  equations 
using  the  Modified  Method  of  Block  Cyclic  Reduction  which  is  described 
in  Appendixes  A  and  B. 
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3.1  SMALL  CHORD  AIRFOIL 

In  the  first  step,  a  simple  vortex  is  chosen  to  represent  the  lift 
model  as 


JL  _  -r  -1  y 

♦m  -  2—  tan  x  (14) 

A  solution  for  the  case  of  a  wall  with  a  uniform  porosity  distribution 
has  been  obtained  to  check  with  the  known  analytical  solution  case  and 
is  shown  in  Fig.  2.  The  second  case,  computed  for  an  inverse  gaus- 
sion  distribution  of  R//3,  is  compared  with  results  obtained  by  the 
approximate  method  (Ref.  2)  in  Fig.  3.  The  agreement  between  the 
results  using  the  proposed  technique  and  previous  solutions  for  the 
above  cases  indicates  that  the  accuracy  of  the  present  numerical  solu¬ 
tion  is  satisfactory. 


-  ANALYTICAL  METHOD 

•  PRESENT  NUMERICAL  METHOD 


Figure  2.  Comparison  of  block  cyclic  reduction  and  analytic  solutions 
for  walls  with  uniform  porosity  distribution. 


3.2  OPTIMUM  POROSITY  DISTRIBUTION 

The  ideal  porosity  distribution  for  a  tunnel  wall  is  defined  as  that 
which  induces  no  lift  interference  anywhere  in  the  test  section.  In  the 
mathematical  sense,  the  upwash  interference,  9 <t>  /9y,  vanishes  every¬ 
where;  or  the  interference  potential  is  a  trivial  solution  of  the  system 
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-  APPROXIMATE  METHOD 

•  PRESENT  NUMERICAL  METHOD 


Figure  3.  Comparison  of  block  cyclic  reduction  and  approximate  solutions 
for  walls  with  inverse  gaussian  porosity  distribution. 


of  Eqs.  (4).  (5).  and  (6).  This  solution  can  be  obtained  by  observation 
as  the  right-hand  side  of  Eq.  (5)  becomes  zero  and  substituting  Eq.  (4) 
then 


R(x)/6 


(15) 


=  x 

However,  the  porosity  parameter  for  the  perforated  wall  R//3  can  only 
have  a  positive  value  because  the  mass  flow  is  always  from  the  high- 
pressure  to  the  low-pressure  side.  Thus,  a  distribution  of  R(x)/j3  is 
selected  and  shown  in  Fig.  4  and  Table  1  denoted  by  Configuration  C  as 


R(x)/6 


1  x/3h 
\  0 


x  >,  o 

x  $  o 


(16) 


to  evaluate  the  interference.  The  lift  interference  factor  for  Configura¬ 
tion  C  has  been  calculated  and  is  shown  in  Fig.  5.  The  interference  fac¬ 
tors  for  three  additional  Configurations  D,  E,  and  F  (Fig.  4  and  Table  1) 
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Figure  4.  Wall  configuration  with  various  porosity  distributions. 


Table  1.  Wall  Porosity  Distribution,  R/0  for  Various  Configurations 


O  00 
0  25 
0  50 
075 
I  00 
I  .25 

1  50 

1.75 
2.00 

2  25 
2.50 

2.75 
300 
3.25 
4.00 
5.00 
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0  250 
0  500 
0  750 
I  000 
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1.500 
I  750 
2.000 
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3250 
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17.000 
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0. 250 
0.500 
1.000 
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1.750 
2.000 

2.250 
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2.750 
3000 
3250 
4.000 
5.000 

17.000 


0.000 


0.250 

0.500 

0.650 

1.200 

1.550 

2000 

2.250 

2.500 

2.750 
3.000 

3.250 
4.000 
5000 

17.000 


0.000 


0  200 
0  4  50 
0.700 
0.975 
I  .300 
1.600 
1 .950 
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Figure  5.  Lift  interference  on  a  small  chord  airfoil  in  tunnels 
with  various  wall  configurations. 


with  a  slight  variation  from  Configuration  C  have  been  calculated  and 
are  presented  in  Fig.  5.  Also  plotted  in  Fig.  5  are  results  for  walls  with 
uniform  (Configuration  A)  and  inverse  gaussian  (Configuration  B)  poros¬ 
ity  distributions.  It  seems  that  Configurations  C  and  D  give,  overall, 
less  interference. 


3.3  FINITE  CHORD  AIRFOIL  CASE 

For  a  finite  chord  airfoil  with  camber  and  incidence,  a  discrete 
distribution  of  vortices  can  be  used  as 


4> 


m 


21 

2tt 


a  .  (£  • ) 
1  J 


(17) 
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The  results  for  the  NACA  64- series  airfoil  with  a  chord  C  =  0.  5/j3h 
are  presented  in  Fig.  6  and  indicate  that  Configurations  D  and  E  exhibit 
the  most  satisfactory  distribution  of  porosity  to  obtain  the  minimum 
interference  factor. 


Figure  6.  Lift  interference  on  a  finite  chord  airfoil  in  tunnels 
with  various  wall  configurations. 


3.4  EFFECT  OF  TEST  SECTION  LENGTH 

Most  analytical  approaches  in  wind  tunnel  theory  have  assumed  the 
length  of  test  section  to  be  infinite  for  mathematical  simplicity.  The 
effect  of  test  section  length  on  the  lift  interference  is  of  interest  since 
the  actual  tunnel  test  section  length  is  usually  about  two  to  three  times 
the  test  section  height.  The  versatility  of  the  present  approach  can  be 
applied  to  examine  the  effect  of  test  section  length.  For  the  uniform 
porosity  distribution  case,  the  comparison  of  lift  interference  of  a 
finite  test  section  as  -2  <  x//3h  <  3  (upstream  and  downstream  regions 
using  solid  walls)  with  the  infinite  test  section  is  shown  in  Fig.  7  and 
indicates  the  effect  on  the  interference  in  the  region  x/ |3h  >  2.  It  can 
be  seen  that  the  assumption  of  an  infinite  length  test  section  for  the 
calculation  of  interference  in  the  neighborhood  of  the  model  appears 
reasonable. 
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INFINITE  LENGTH  TUNNEL 

FINITE  LENGTH  TUNNEL  (-2  5  x/£h  <3) 


Figure  7.  Effect  of  test  section  length  on  lift  interference. 


4.0  CONCLUDING  REMARKS 


An  efficient  numerical  scheme  has  been  developed  by  a  modification 
to  the  Block  Cyclic  Reduction  Method  for  computing  lift  interference  in 
a  wind  tunnel  with  an  arbitrary  distribution  of  wall  porosity.  A  compar¬ 
ison  with  other  available  analytical  and  approximate  solutions  has  demon¬ 
strated  the  accuracy  of  the  present  numerical  method.  The  optimum 
porosity  distribution  to  minimize  interference  is  obtained  by  the  varia¬ 
tion  of  the  ideal  mathematical  configuration  which  produces  exact 
interference-free  condition.  The  minimization  of  interference  is  also 
presented  for  a  finite  chord  airfoil  in  the  optimum  wall  configurations. 

The  optimum  wall  porosity  configurations  have  been  calculated  for 
both  a  simplified  ideal  airfoil  and  a  finite  chord  airfoil.  The  effect  of 
test  section  length  has  been  also  studied. 
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APPENDIX  A 

METHOD  OF  BLOCK  CYCLIC  REDUCTION 


Consider  the  problem  of  solving  the  finite  difference  analog  to 
Laplace's  equation  with  the  boundary  conditions 


<(>  (-x*,y)  =  <f>(x*,  y)  =  0  -1<  y  <  1 


<J>  (x,  +1)  =  g~  (x) 


-X*<  X  <  X* 


and  where  g  (x)  are  given  functions  with 

g1 (x*)  =  g' (-x*)  =  0 . 


(A-l) 


It  is  well  known  that  replacing  Laplace's  equation  on  the  region 
R  by  a  centered  second  difference  approximation  and  imposing  the 
boundary  conditions  given  in  Eq.  (A-l)  yields  the  problem 


D  $=  y 


(A-2) 


where  D  is  the  (2M-1)  (N-l)  x  (2M-1)  (N-l)  real  symmetric  matrix 
which  has  the  block  tridiagonal  form 


D  -  (I,  A,  i)2m-1  x  2M-1 

and  A  is  the  matrix  defined  in  Eq.  (9). 

'  The  vector  q>  will  be  given  in  partitioned  form  as 

$  =  co1(^m-1'  ?M-2'  ?1-m)  ’ 

Likewise,  the  vector  y  is  given  by 

y  =  COl  (“  r  0  ,  .  ...  t  0  t  — 

In  their  description  of  Block  Cyclic  Reduction,  Buzbee,  Golub,  and 
Nielson  (Ref.  6)  first  write  Eq.  (A-2)  as 
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A!m-1  +  +M-2 

~M-1 

(A- 3a) 

A^j  +  =  o 

j  =  0#+l#...r+ (M—  2 ) 

(A- 3b) 

*2-M  +  A?l-M 

=  ~1-M  ‘ 

(A- 3c) 

Then  for 
written 


=  i- 1,  SL,  j2  + 1  where  1!  =  -M+2,  .  .  . ,  M-2,  Eq.  (A-3b)  can  be 


h+2  +  At*+i  +  t* 


it+i  +  Ah  +  tt-i 

h  +  A!n-i  +*ji-2 


0 

0 

0  . 


Then  multiplying  the  middle  equation  by  -A  and  adding  the  three 
equations  yields 


$1+2  +  (2l"A2)h  +  *£-2  =  0 

for  £  =  -  M+2f  ...  ,  M-2  . 


Buneman,  as  described  by  Hockney  (Ref.  7)  proceeds  by  these  steps 
and  then  reapplies  the  method.  Thus  with 

=  2 I -A2  , 

$SL+ 4  +  Al$Z+2  +  $1  =  0 

$1+2  +  A1 $Z  +  $£_2  ’  =  ° 

$1  +  Al$£-2  +  $1-4  =  0 


where  £  =  -  M+4 ,  ...  ,  M-4 

and  again,  multiplying  the  middle  equation  by  -A,  and  adding  yields 

$£+4  +  {2I_A1)$£  +  h-4  =  0  * 

Then  repeating  the  process  of  cyclic  reduction  recursively,  Buneman 
obtains  for  the  i**1  recursion 
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<t>  j  +  A.  0  .  +  0  •  =  0 

~  j+21  x~3  "*  j-21 


A.  =  21-Af  , 

x  1-1  ( A-4) 

A0  -  A 

Hence,  when  j  =  0  and  i  =  k  there  results 

$M  +  Mo  +  $-M  =  0 

so  that 

*0  “  '  ^'tM  +  !-M>  • 

djyj  and  0_j^  are  known  values  from  Eq.  (A-l);  hence,  0 g  may  be 

found  by  inverting  an  (N-l)  x  (N-l)  matrix.  Once  0g  is  known,  the 
method  may  be  repeated  on  the  regions 


«u- 


(x,y)  -  x*  <  x  <  xi 


0  <  y  <  1 


and  R_  = 
L» 


(x#y)  -  x*  <_  x  <_  x*j 


-  1  <  y  <  0 


solving  for  0M  and  0 

~2  ~T 


These  steps  are  repeated  until  all  the  vectors  0^-2  =0,  ±1,  . . .  ,  ±(M-1) 
are  found.  Each  step  requires  finding  the  solution  to  N-l  linear 
equations. 
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APPENDIX  B 

MODIFICATION  OF  BLOCK  CYCLIC  REDUCTION 


The  set  of  finite  difference  equations  (Eqs.  (8),  (11),  and  (12))  was 
developed  in  Section  2.  2  and  is  given  by 

d> .  , ,  +  Ad> .  +  sb  .  i  =  0 
~J+1  ~J-1 

j  =  0,  +  1,  ...  #  (M-l)  (B- 1) 

t*m-i  “  e*m  +  !+  (B'2a) 

t»i-m  =  e*-m  +  f  •  <B'2b) 


It  will  be  shown  that  the  vectors  and  can  found  as  the  solu¬ 
tion  to  a  system  of  2(N-1)  linear  equations. 

At  this  point,  a  change  of  notation  will  be  made  for  convenience. 

Let 

Ya+M  =  h  *  -  -M,  ...  ,  H  .  (B-3) 


Applying  Eq.  (B-3)  to  Eq.  (A-4)  results  in 


V  .  +  A.  V. +  V  .  =  0 

~  j+21+M  1  M  '  j-21+M 


(B-4) 


Theorem 


Yjt+1  =  Fn  Yz  +  Gn  Y£+2n 


(B-5) 


where 


n  =  1,  2 1  ...  r  K+ 1 


and  £  is  such  that 


%  >  0  and  2n  +  l  <  2M  , 
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where 


(B-6) 

(B-7) 

and  where 


F„  -  Fn-1  *  V-l  Vl  n=  2>  •••  '  K+1 


G_  =  -  G_  ,  A_.  n  =  2,  ...  ,  K+1 
n  n-x  n— x 


Proof 


In  Eq.  (B-4)  let  j  =  1  -  m  +J2  and  i  =  0. 


Then 


Y2+1  +  Ao  Yi+i  +  Y*  -  0 


hence 


Ywi  -  -  Ao1  <Yi  +  Yi+2> 
-  Fi  Yi  +  Gi  Y1+2 


so  the  theorem  holds  for  n  =  1.  These  steps  would  complete  the  proof 
for  K  =  0  so  now  assume  K  >  0  and  suppose  Eq.  (B-5)  holds  for  n  =  L, 
L  =  1 . K. 


Then 


Y*+i  ■  fl  Y»  ♦  =L  YU2l 


(B-8) 


where  S.  is  such  that 

£  >  0  and  2L  +  £  <  M  . 

In  Eq.  (B-4)  then,  let  j  =  2^  -  2^  +  H  and  let  i  =  L. 

Then 

Y„l+i  .  +  v  Y,l  .  +  Yi  -  0 
2  +£  2+1 
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or 


v+. 


-1 


■  -  *1  (*» +  VO  • 

Substituting  Eq.  (B-9)  into  the  inductive  hypothesis  Eq.  (B- 

Yji+i  -  lL  ll 

=  (fl  -  gl  *£l)  Y*  - 


v-  =  p' v"+gL[-^  +  VO 


Gl  V  Y2l+i+1 


fl+i  Yi  +  gl+i  Y2l+i+a 


and  hence  the  proof  is  complete. 

Then  with  n  =  K+l  and  Si  =  0.  Eq.  (B-5)  becomes 


Yi  “  fk+i  Yo  +  gk+i  Y2k+i 
=  fk+i  Yo  +  gk+i  Y2M  * 


In  a  similar  manner,  it  may  be  shown  that 


Y2M-1  fk+i  Y2M  +  gk+i  Yo  * 


Applying  Eq.  (B-3)  to  Eqs.  (B-10)  and  (B-ll)  gives 

Jl-M  =  FK+1  $-M  +  GK+1  $M 
$M-1  =  FK+1  $M  +  GK+1  ~-M  ' 


Substituting  Eqs.  (B-10)  and  (B-ll)  into  Eqs.  (B-2a)  and  (B-2b) 
tively,  yields 


T[FK+1  +  GK+1  $-m]  E$M  +  ~ 
T[FK+1  ~-M  +  GK+1  $m]  =  E$-M  +  ~ 


(B-9) 
8)  gives 


(B-10) 

(B-ll) 

(B-  12) 
(B-  13) 

respec- 

(B-  14) 

(B-  15) 


24 


AEDC-TR-75-98 


which  may  be  written  in  Block  Matrix  form  as 


TFK+1  "  E  gk+i 

1 

»-©- 

s 

_ 1 

1 - 

+ 

•H, 

gk+i  tfk+i  "  E 

$-M 

\3'1 

It  will  be  noted  that  Eq.  (B-16)  is  a  linear  system  of  2(N-1)  equa¬ 
tions  which  can  be  solved  for  and  Once  these  vectors  are 

known,  the  problem  becomes  one  of  the  Dirchlet  type  which  can  be 
solved  by  the  methods  of  Appendix  A, 

In  Theorem  I,  let  n  =  K  and  SL  ~  M;  then  Eq.  (B-5)  becomes 

Y  M+l  =  FK  Ym  +  GK  Y2M  *  (B- 17) 

Also,  writing  Eq.  (B-4)  with  j  *  0,  i  =  K  and  again- noting  that  2K'M 
results  in 


or 


Y2m  +  ak  Ym  + 


Yo 


o 


Ym  "  -  ak1(Yo  + 


Y2M> 


(B-  18) 


Substituting  Eq.  {B-18)  into  Eq.  (B-17)  gives 

Ym+1  =  Fk[“  AK  (Yo  +  Y2m)]  +  GK  Y2M 
=  (gk  fk  *K  )Y2M  -  fk  ak  vo 

Letting 

S  =  Gr  -  Fr  A^1  and  W  =  -  FR  A^1 

gives  Ym+i  -  s  Y2m  +  w  Yo 

and  then  by  use  of  Eq.  (B-3) 

4l  “  S*M  +  W*-M  ' 
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In  a  similar  manner  it  is  found  that 

+-i  =  w!m  +  s!-m  • 


Then 


h  "  t- 1  =  (s_w)  (^M  "  *  (B-19) 

Subtracting  Eq.  (B-15)  from  Eq.  (B-14)  results  in 

(tfk+i  -  tgk+i  -  e>  %  -  ♦.«>  -  <?+  -  r> 


so  that 

♦  l  "  $_!  =  <S-W)  (TpK+i  -  TGk+1  -  E)'1  (f+  -  f“)  (B- 20) 

hence,  0^  -  can  be  computed  by  solving  an  N-l  system  of  equa¬ 
tions  and  the  interference  factor  in  Eq.  (13)  is  obtained. 
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APPENDIX  C 

METHOD  OF  EVALUATION 


The  evaluation  of  the  lift  interference  by  use  of  Eq.  (B-20)  is 
greatly  hindered  by  the  number  of  operations  required  to  evaluate  the 
recursion  matrices  Fk+j  and  .  However,  these  computations 

may  be  greatly  simplified  by  a  simple  application  of  induction  and  it 
is  shown  by  use  of  Eqs.  (B-6)  and  (B-7)  that 


and 


so  that 


Gk  =  (-1)K(Ak_1  ...  .  Ax  AQ)  1 


K 


Fr  =  I  g5 


1=1 


K+l 


FK+1  “  Gk+1  21 


K+l  GJl  “  GK+1  FK  ’ 


Since  S  -  W  =  Gj^ ,  Eq.  (B-20)  becomes 

h  "  t-i  =  gk(tpk  "  E>”1(?+  ~  !_)  * 


(C-l) 


Now  consider  the  matrix  A  given  by  Eq.  (9). 


Define  the  matrix 

D  =  [°'  di'  °]N-1  x  N-l 


where  d^ 

A 

Define  A 


where 

A 

and  c 

l 


1  and  dj+1  =  aj  d./Cj+1,  j  =  1,  . ..,  N-2.  (C-2) 

A  [c.f  b±, 

b±  i  =  1, 

ai-i  =  V°i  ai-i 


ai]  N-l  X  N-l 
. . . ,  N-l 


(C-3) 
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A 

Then  A  is  a  real  symmetric  matrix  for  which  many  well  known  computer 
programs  can  be  used  to  compute  the  eigenvalues  and  eigenvectors  X  j 
and  xj  .  j  .1,  ,  N-l. 

A 

From  Eq.  (C-3),  A  and  A  are  seen  to  be  similar  matrices  hence 
Xj  and  D"l/2Xj  form  an  eigenvalue,  eigenvector  pair  for  A. 

Let  X  be  the  unitaryAmatrix  whose  columns  are  the  set  of  ortho¬ 
normal  eigenvectors  of  A. 

Then  X-1  A  X  =  A 
where  A  -  [o.  o]  ^  x  „.1  . 

Hence  Q"  *  A  Q  =  A  where  Q  =  D~  ^  ^  X. 


Now  suppose  there  exists  a  matrix  which  diagonalizes  A^  so  that 
Qk1  Ak  Qk  =  Ak  • 


Then  since 


=  21 


Ok1  flk+l  Qk  “  Ok1  21  Qk  -  Ok1  Ak  Qk 


Hence,  Qjj  diagonalizes  both  and  Qk+i.  Then  dropping  the  subscript 
on  Q  gives 

Q- 1  Ak  0  -  Ak 

so  that  from  Eq.  (C-l) 

Gk  =  (-l)k  (Ak_1 - Ao)_1 

=  (-Dk  (Q\_x  Q”1  Q\_2  Q"1  ...QAqQ"1) 

=  (-l)k  (Q\_! - Ax  A0  Q"1) 
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But  Q"1  =(D"1^2X)"1  =X_1  D1/2 
and  since  X  is  a  unitary  matrix, 

Q-1  =  XT  D*  .  So  finally, 

Gk  =  (-1) k  {D-^  X  Xk-1  - X±  Aq  XT  D^)  . 

Hence  the  work  of  computing  G^  is  simplified  since  most  matrices 
involved  are  diagonal  matrices. 

It  will  be  noted  that  and  are  functions  only  of  the  mesh  and 
hence  the  interference  may  be  computed  via  Eq.  (C-2)  for  many  dif¬ 
ferent  porosity  distributions  T  without  having  to  recompute  these 
matrices. 
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APPENDIX  D 
COMPUTER  PROGRAM 

Program  Description 


MAIN 

The  main  program  is  for  control  purposes.  SETTUP  should  be 
called  immediately.  TFIX  is  called  whenever  a  new  porosity  distribu¬ 
tion  is  desired.  EVALU8  is  called  to  compute  the  lift  interference. 

In  the  sample  listing,  the  interference  is  computed  for  the  four  poros¬ 
ity  distributions  given  in  Table  1. 

B 


Function  B  is  a  user- supplied  routine  and  is  used  to  compute 


+ 

BT  =  - 

l 


R(V  a»m<xl, 

8  dx 


.  <Xi, 

3y 


This  equation  is  similar  to  the  one  following  Eq.  (9).  It  need  be 
rewritten  only  when  the  model  velocity  potential  is  changed. 


CHOLES 

Subroutine  CHOLES  solves  the  linear  system  given  by  Eq.  (C-2) 
by  the  method  of  matrix  factorization. 


DFFIX 

Subroutine  DFFIX  computes  the  vector  f +  -  f". 

EVALU8 

This  subroutine  evaluates  the  finite  difference  coefficients,  con¬ 
structs  the  matrix  (TF^  -  E)  G^1,  and  calls  subroutine  CHOLES  to 
solve  Eq.  (C-2). 


FFIX 


This  subroutine  constructs  the  matrices  Fk  and  Gk  by  the  methods 
described  in  Appendix  C. 
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HDIAG 

Subroutine  HDIAG  computes  the  eigenvectors  and  eigenvalues  of 
the  matrix  A  defined  by  Eq.  (C-3). 

MESH 

This  user- supplied  subroutine  is  used  to  fill  the  X  and  DX  arrays. 
Note  that  a  value  is  assigned  to  X(o),  namely  X{o)  =  -X*. 


MULT  and  MULT2 

These  routines  perform  the  FORTRAN  matrix  replacements 
B  =  AB  and  A  =  AB,  respectively. 


SETTUP 

This  is  a  user- supplied  routine  used  to  initialize  all  program 
constants. 

TFIX 


TFIX  loads  the  vector  T(x^)  by  calling  TFUNC. 

TFUNC 

This  is  a  user-supplied  routine  used  to  evaluate  the  function. 

TFUNC  (X)  =  R(x)/6 

It  should  be  noted  here  that  in  the  program  the  array  T  contains  values 
of  R//3  and  not  j3/R  as  given  following  Eq.  (9). 
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MAIN 


IMPLICIT  REAL*8(A-H,0-2) 

COMMON  T(49),VEC(49),0UM(49,50) 

COMMON/XX/XDf  XI 50) 

COMMON /I 1/N,M,N1 ,Ml,K 

COMMQN/CONTRL/IC 

CALL  SETTUP 

00  50  ICC- It  4 

IC=  5-  ICC 

CALL  TF1X 

CALL  EVALU8 

00  39  I  =  l»Nl 

39  WRITE  (6,  10)1  ,X(  1 1  «T|  I  I  ,DUM{  I  ,N) 

10  FORMAT ( 5X, • I* ' , 1 3, 5X, • X= • ,F10.3 ,5X , *8ETA/R0  =  « ,016,8,  5X, 
**DELTA=',D16«8) 

50  CONTINUE 
STOP 
ENO 


B 


DOUBLE  PRECISION  FUNCTION  B(J) 
IMPLICIT  REAL*8(A-H,0-Z) 

COMMON /Rl/ 0X0, D  XI  49)  ,0UM(5)  ,XINF,P12 
COMMON /XX/ XD  ,E  X( 50) 

COMMON  T ( 49) 

COMMON/BROWN/GAM 
1= I ABSi J  ) 

X=  EX(  I ) 

Y=  1/ J 

B=(-Y*Tt  I)*Y*X)/(X*X*l.DO)  *GAM/PI2 

RETURN 

ENO 
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C  HOLES 


SUBROUTINE  CHOLES(A,N 
REALMS  A {  101,1021. SUM 

,NV, 101,1 02 , MAT SYM) 
.TEMP 

M=N*NV 

NARD=N  +  1 

IF ( At  1, 11. NE. 0.01  GO 
DO  37  J=2,  N 

TO 

47 

IFIAC  J,  D.EO.O.O)  GO 
IFLIP*J 

TO 

37 

GO  TO  27 

« 

37 

CONTINUE 

GU  TO  54321 

27 

DO  57  K*l.M 

TEMP® A ( 1 FL  IP  ,K  ) 

At IFL IP.K )=A( 1. K) 

At  i,K  )=TEMP 

57 

CONT  INUF 

47 

DO  2  J=2»M 

At  1. J »  =  A( l.JJ/AI 1, 11 

2 

CONTINUE 

DO  6  1=2, N 

DO  7  J=2,M 

IF (MAT  SYM .EQ . 0 IGO  TO 

49 

IF( I— J  J69, 68,67 

49 

IF ( J  «UT. 1 ) GO  TO  69 

68 

K= J- 1 

SUM=0.0 

DO  3  1R=  1, K 

SUM=SUM*A( I, IR)*At 1R, 

J) 

3 

CONT INUE 

At  I,  J  )=AI  I  ,J  I-SUM 

GO  TO  7 

69 

K=  I-  l 

SUM=0.0 

DO  4  IR-l.K 

SUM=  SUM+AC  I,IR)*A(  IR, 

J) 

4 

CONTINUE 

TfTA'I  I.ITTEQ.O.OI  GO  TO  54321 
At  I, J )  =  ( A1 I,J}-SUM)/A< 1,1) 

GO  TO  7 


67  A(  1,J)=AU,  N»At  J«  J» _ 

7  CONTINUE 

6  CONTINUE  _ 

00  52"  NPR6B"=NARD  »M 

DO  52  K=2,N  _ 

I  =  N«-l-K 

SUM«=  0.0 _ 

Ll*I+l 

DO  51  IR*LL  »N  _ 

SUM=SUM«-"A<  I,IR)*A< IR.NPROBI 

51  CONTINUE 

At  I,NPROB')*AlI  ♦NPROBI-SUM 

52  CONTINUE 

GO  TO  12345 

5432 1  _N«- 1  _ 

12 345  RETURN 
END 
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OFF  l  X 


SUBROUTINE  UFFIX 
IMPLICIT  K  EAL*8l A— H,0-Z ) 
CGMMON  T (  49)  ,  VEC l  49) 
COMMON / 1  l/N,M, N1 
COMMON /R  1/0X0,0X149)  ,OY 
DU  l  I  =  l ,N  l 

I  VEC(i)=-DY*(B( i )-8 I— I ) ) 
RETURN 
END 


E VALUd 


SUBROUTINE  EVALU8 
IMPLICIT  RfcAL*8<  A-H,0-Z) 

COMMON  T( 49 ) ,VEC{ 49) ,0UM(49,50I ,FKI49,49)  , APROD 149 , 49 1 
COMMON /R  l/DX0,DX(49)  , 0 Y  ,B ETA  , PO ,  fc,CS  » XINF ,PI2 ,Y ,DY2 
COMMON /I l/NfM,NI , M 1 » K , I  0 1 «1D2,IZ,N2 
COMMON /ABC /A ( 49 ) ,6(49) ,C(49) 

100  CALL  OFF  IX 
00  l  1  =  1,  Ni 
DO  1  J=1*N1 
1  DUMC  I,  J)=FKl  I,J) 

Q  A=DY2*l  DX<  D-DXI  I  Z  ) )  / (  DX(  1 )  *DX (  I  Zl ) 

Pl=DY2*DX<  1Z)/(DX( l)*(UX(  1UDXCIZ)  )) 

DUMl  l,  i)=DUMC  l,  1)-.5D0*(  TCI)  *01-6(1)  I 
OUMI  1,  2)=D'JM<  l,2)-.5D0*(  TCI) *P1-A(1)  ) 

DO  65  1=2, N2 

QI=DY2*(0X(  I  )-DX(  I-i)  1/(0X111*0X11-1)  ) 

PI  =  DY2*DX(  I—  1 )  /  (  DX  (  I)*<DX(n+DX(I-ll  )) 

Ri=-DY2*DX( I)/(DX< I-1J«(DX|1 )+DX(l-l*)i 
DUM( I, I-1)=DUMC I 5D0*(T( I ) +RI-CII) ) 

DUM(  I,  I  )=OUM  (  I,  I  )-.500*<T(l)*QI-B(I)  I 
65  DUMC  I,  I*  1 1  =DIJMI  i ,  l  ♦  L )-.  5D0*(  T  C 1 )  *P1-A(  1 )  1 
Rl=-DY2*DX(N1)/(DXC  N2 )  *  (DX  C  Nl  J+DXCN2 )  ) ) 

0  l=DY2*(  DX(Nl)-DX(  N2) ) /(DXCN1) *UX(N2) ) 

DUMC  N 1 ,N2) =DUM( N l,N2)-«5D0*( TC Nl I *Rl -C CNl ) ) 

DUMC  Nl, Nl) =DUM( N1,N1)-.5D0*( T( Nl ) *Ol-B( Nl ) ) 

CALL  MJLT2CDUM, APROD) 

DO  20  1=1, Nl 
20  DUMC I,N)=VECC I ) 

CALL  CH0LESCDUM,N1,1, 101,1 D2,0) 

DO  4  I  =  1  ,N  l 

4  DUMC I,N)=DUMC I ,N ) / ( 2. DO*DY) 

RETURN 

END 
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FF  IX 


SUBROUTINE  FFIX 
IMPLICIT  REAL*8(  A-H,0-Z) 

COMMON  T(49),VEC(49),TMP<49,50) ,FKl49,49) ,APR0D(49 , 491 
COMMON /Rl/ OXO, OX( 49) ,0 Y^DUMI 5) ,P12 
DIMENSION  D( 49) ,0GI49) ,0F(49) 

COMMON /I l/N,M,Nl,Ml,K 
COMMON /ABC /A (49) ,B(49) ,C(49) 

DY2=2.D0*DY*DY 
DO  20  1=1, N1 

A(I)  =  DY2/CDXCI  )*IDX(  I  J  +  DXlI-l)  )  ) 

B 1 1 )=-2.D0-DY2/<DX{ I ) *0X1 1-1) ) 

20  CU)  =  DY2/(DXU-U*{0XU  )+DX(I~l)  )) 

0(  L)=1.00 
N2=N  I- 1 
DO  30  1=1, N2 

30  0(  I«-1)=A1  I  1*0 (  I  l/CU  +  l) 

00  60  1=1, N1 
DO  60  J= 1 *N  1 
60  FK  I  I ,  J  )=  0.00 
FK( 1, l)=B( 1) 

DO  40  1  =  2, N1 
FK< I, I  )=B (  I  ) 

FK{  I,  I- 1  )  =  DSQRT (C  (  I  >*A(  I-D) 

40  FKI 1-1,1 )=FK< 1,1-1) 

I EGN=0 

CALL  HDIAG(FK»N1»IEGN,TMP,NRN) 

DO  50  1=1, N1 

VEC(  I  )  =  1  .DQ/DSQRT  (0(D) 

50  D(  I  )=FK(  1,1) 

DO  1  1=1, N1 
OF  I  I )=  0.00 

1  DG(  I  )=  1.00 
00  2  IK= l,K 
DO  2  1=1, N1 

DG(  I  )=DC  I  )*OG(I  ) 

DF( I )=I- 1) ** IK/OGI I > ♦OF  III 

2  DI I  )  =  2  .DO-DI  I )  **2 
EE= (  —  1  )**K 

DO  10  1=1, Nl 
10  DG(  I  )=EE*DG(  I  ) 

DO  501  J=1,N1 
00  501  1=1, Nl 
APRODI  I,  J)=TMP(  I,J  )*DG(  J) 

501  FKI  I,  J)=TMP(  I,J)*DF(  J) 

DO  502  1=1, Nl 
DO  502  J= 1 , N 1 

APRDOI i, J)=VEC( I ) * APROD ( I , J) 
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502  FK(  I,J  l=VEC<  I>*FM  [  ,  J) 

DQ  503  1  =  2, N1 

IM  1=1-1 

DO  503  J=i,IMl 
TEMP=TMP<  I, J) 

TMPi I, Ji=TMP( J, 1) 

503  TMPC J,  I  )=TEMP 

CALL  MULT2(APR0D,TMP) 

CALL  MULT2(FK,TMP) 

DU  504  J= l ,N l 
TQM= l • DO/VEC ( J ) 

DO  504  1  =  1, Nl 

APRO  D( I , J )=APRQD( I ,J)*TOM 

504  FK  (  i,  J  }=FK  (  I  ,  J  )  +  TOM 
RETURN 

END 
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HD  I  AG 


SUBROUTINE  HD I  AG  (H,N,IEGEN,U,NR> 

C  SUBROUTINE  hDIAG. _ 

C 

C  PROGRAMED  BY  F.  J.  CARBATO  AND  H.MERW1N  OF  THE  HIT 
C  COMPUTAflUN  CENTER^  '  1 

C _ 

C  THIS  SUBROUTINE  COMPUTES  THE  EIGENVALUES  ANO  EIGENVECTORS 

C  OF  A  REAL  SYMMETRIC  MATRIX,  H,  OF  ORDER  N  I  WHERE  N  MUST  BE  LESS 

C  THAN  5l»,  AND  PLACES  THE  EIGENVALUES  IN  THE  DIAGONAL  ELEMENTS  OF 

C  THE  MATRIX  H,  AND  PLACES  THE  E1GENVE CTGRS  (NORMALIZED  >  IN  THE 

C  COLUMNS  OF  THE  MATRIX  U.  1EGEN  IS  SET  AS  1  IF  ONLY  EIGENVALUES 

C _ ARE  DESIRED, AND  IS  SET  TO  0  WHEN  VECTORS  ARE  REQUIRED.  NR  CUN- 

C  TAINS  THE  NUMBER  OF  ROTATIONS  DONE. 

C _ _ 

C  H,  N,  IENGEN,  U,  AND  NR  OF  THE  ARGUMENT  LIST  ARE  DUMMY  VARIABLES 

C  ANO  HAY  BE  NAMED  DIFFERENTLY  IN  THE  CALLING  OF  THE  SUBROUTINE. 

C 

C _ SUBROUTINE  PLACES  COMPUTER  J  N  THE  FLQAT I  N«  TRAP_MGDE _ 

C  THE  SUBROUTINE  OPERATES  ONLY  ON  THE  ELEMENTS  OF  H  THAT  ARE  TO  THE 

C _ RIGHT  OF  THE  MAIN  DIAGONAL.  THUS,  ONLY  A  TRIANGULAR _ 

C  SECTION  NEED  BE'  STORED  IN  THE  ARRAY  H. 

_ IMPLICIT  REAL*8( A -H,Q-Z> _ 

DIMENS  ION  HI  49, 491 ,U! 49,491 , X (49) ,IQt49> 

_  2 _ FORMAT ( 14H  MAX  OFF  DlAb=,F14.7,3hNR»,I3) 

2001  FORMAT  I  LX, BE  15. 8) 

2002  FORMAT  I  IBH  ORTHOGONAL  MATPJ_XJ_ _ 

2003  FORMAT! 15H  ROTATED  MAtRIX) 

_ SIGN! X  l,X2)=0SI GNC  XI , X2| 

SURTIXI-DSQRTI X) 

ABSI XI=DAB SI  XI 

IF(  IEGEN.NE.O)  GO  TO  15 

10  00  14  1=1, N 
00  14  J=1 ,N 

_ IF!  I-J.NE.O)  GO  TO  12 _ 

11  U(  I, J )-l .0 
GO  TO  14 

12  U(1,J)=0.0 

14  CONTINUE 

15  NR  *  0" 

_ IFCN-  1  .LE.01  GO  TO  1000 _ _ 

C  SCAN  FOR  LARGEST  OFF-DIAGONAL  ELEMENT  IN  EACH  ROW 

C  X(I»  CONTAINS  LARGEST  ELEMENT  IN  I TH  ROW 

C  IQ!  1)  HOLDS  SECOND  SUBSCRIPT  OEFlNlNu  POSITION  OF  ELEMENT 

17  NMI1=N-1 _ 

DO" 30'  I-l.NMI  I 

_ HI)  »  0.0 _ 

IPl  1=1+1 
DO  30  IPL  1,N 

1-ABSI HI  I  ,J>  ).GT.  0.0)'  GO  TO  30 
20  XlI)=ABS(HII,jn 

mcnvj- 

30  CONTINUE _ 

C  SET  INDICATOR  FUR  SHUT-OFF. RAP=2**-2 7 ,NR»NO.  GF  ROTATIONS 
RAP* 7.4505805960-9 
HDTE ST* 1.0038 

C  FIND  MAXIMUM  OF  XII)  S  FOR  PIVOT  ELEMENT  AND 
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C  TEST  FOR  END  OF  PROBLEM 
40  DO  70  1=  It NMI I 

iFli-l.LE.OI  GO  "TO  60 

_ I_FI_XM_A_X-X I  I ). JE.O.O)  uU  TO  70 

60  XMAX=XUI 

_ IP  IV= I_  _ 

JP  IV= 1QI  I  ) 

70  CONTINUE 

C  IS  MAX.  XII)  EQUAL  TO  ZERO*  IF  LESS  THAN  HOTEST,  REVISE  HDTEST 

_ 1FC  XHAX.LE.O.OI  GO  TO  1.002 _ 

80  IF(HOTEST.LE.O.O)  GO  TO  90 

_ 8  5 _ I F  {  XMAX-HDTE  SJ[  .GT  .0.0)  GO_TO_ 148 _ 

90  HDIHIN  =  ABS<  H(  It'll  '  f . 

_ DO  110  1  =  2. N  _ _ 

IF  (HD  MIN-  A8S(  H(  f,  I  )  l.LE.  0.0)  GO  TO  110 

100  HDIM jN  =  ABS_|  H( 1 ♦ I )  ) _ 

110  CONTINUE 

_ HOTEST  =  HDIM1 N*RAP _ _ _ 

C  RETURN  IF  MAX.H( I , JILfcSS  THANI2* *-27 ) ABS I  HI K , K l-Ml N I 
IF ( HOTEST- XMA X.bE .0. 0)u0  TO  1002 
148  NR  =NR*1 

C  COMPUTE  TANGENT,  SINE  AND  COSINE .H <1.1* ■HlJ, J ) _ 

150  T  ANl»=  S  luN  I  2.0,  (H|  (P1V,  1PI  V)-H(  JPI  V,  JPI VII) *H 1 1 P I V.  JPIV)/  (ABS  IHC I 
IP  IV,  IPIVI-HI  JPIV,  JPIV)  )  ♦StiRTI  (HIIPIV,IPIV)-H(  JPIV,  JPIV)  )  **2*4.0*  H 

21  I  P  I V  i  JP  I V  j  **2 )  )  .  -  -  -  - 

_ COS  INF  =  l  .O/SQR  T (  l.  0*  TANG**2_I _ 

SINE=TANG*COSINE 

_ HI  1=  HI 1PIV, IP  1 V  ) _ 

H(  IP  IV,  IP  1  V)=COSlNE**2*(HlI*TANta*(2.  *H(I  PI V,  JPIV) ♦TANG*M JPIV, 

_  1JP IV)  )  I  _  _ _ _  _ 

H<  JP I V. JPI V1=C0$INE**2*(H ( JPI V. JPI VI -TANG*(2. *H( I PIV,JPIV)-TANG*H 

_ 1III) _ 

"hi  iPiv.jp iv )=6.o 

C  PSEUDO  RANK  THE  EIGENVALUES _ 

C  ADJUST  SINE  AND  COS  FOR  COMPUTATION  OF  H(IK)  ANO  U(1K1 

_ _ I F ( H(  IPIV,  IPIV)-H( JPIV, JPIVI.GE. 0.0)  GO  TO  153 

152  HTEMP  =  H( IPIV, IPIV) 

_ h(  IPIV,  IPIV)_=_H(JPI  V,JPI_V) _ 

HI  JPIV. JPIV)  =  HTEMP 

_C _ RJ COMPUTF.  SINE  AND  COS _ 

HTEMP  =  SIGN! 1.0,  -SINE)  *  COSINE 
COSINE  =  ABS  (SINE) 

S  INE HTEMP" 

153  CONTINUE 

C  INSPECT  THE  IQS  BETWEEN  1*1  AND  N-l  TO  DETERMINE 

C  WHETHER  A  NEW  MAXIMUM  VALUE  SHOULD  BE  COMPUTED  SINCE 

C  THE  PRESENT  MAXIMUM  IS  IN  THE  I  OR  J  ROW. 

DO  350  I  =  1.NM11 

IFI  I-IPIV.EQ.OI  GO  TO  350 

_ IF( I-IPIV.LT._0  )  GO  TO  210 _ 

200  IFI  I-JP1V.EQ.  O'  I  toO  tq  350 

210  IFI  IQ  I  I  )-  I P  I  V.E  Q.  0)  GO  TO  240 
230  IFI  IQ( II-JPIV.NE.  0  I  GO  TO  350 
240  K  =  IQ  (  1  IS 
250  HTEMP  =  HC'I.K) 

HI  I ,K  )  =  0.0 
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IPL1  =  1+1 

_ X(I)  =  0.0 _ _ 

C  SEARCH  IN  DEPLETED  ROM  FOR  NEM  MAXIMUM 

_ DO  320  J  =  1PL1.N _ 

l  F  {  X  (  I )"-  A  8  S  (  HII.JI  ).GT.  0.0)  GO  TO  320 

300  XC  I)  *  A  0  S  ( till  .  J)  ) _ 

10  If)  =  J 

320  CONTINUE _ 

HII.K)  =  HTEMP 

350  CONTINUE _ 

X<  IP  I  V  )  =  0.0 

_ X  { JP  I V  )  =  0  .0 _ _ 

C  CHANGE  THE  OTHER  ELEMENTS  OF  H 

_ 00  530  I  =  1  *  N _ 

1 F I  I—  IP  I V.EQ •  0  )  GO  TO  530 

_ 1F1 I-1P1V.GT.  0  )  GO  TO  420 _ 

370  HTEMP  =  HIl.IPIV) 

_ HI  I.JP  IV  )  =  COSINE *HTEHP  »  SlNE_*rt|I,  JPIV)  _ 

IF  1  XU)  -  ABSI  HIl.IPIV)  ) . GE.  0.6  )  GO  TO  390 

380  XI  1)  _■ _ ABSI  HI  1 ,  IP  IV)  ) _ 

10 1 i)  =  IP  I  V 

390  HI ItJPlV)  »  -SINE+HTEMP  ♦  COSINE *HII .JPIV) _ 

IF  I  Xtn  -  ABSI  HU. JPIV)  ).GE.  0.0  )  GO  TO  530 
400  XI  I)  =  ABSI  HI  I, JPIV)  ) 

IQ(I)  =  J> IV  " 

GO  TO  530 

420  Iff  I- JP I V.EQ.  0  )  GO  TO  530 

_ IF  I  1— JP I  V. GT»  O  )  GO  TU  480 _ 

430  HTEMP  S  HI  IP  IV. I ) 

HI  IP  I V  •  I  )  =  COS  I NE*HTEMP  ♦  SINE*HI 1 , JPIV) 

IF  I  XI  IPlVr-  ABSI  Hi  IPIVtl)  )  .GE.  0.6  )  GO  TO  450 
440  XIIPIV)  =  ABSI  H(IPIV.I)  ) 
fQ!  IP  TV  )  =  I 

450  HI1.JPIV)  =  -SlNE*rlTEMP  ♦  COSI NE  »H  II  .  JP1  V) _ 

IF{  XII)  -  ABSI  HU. JPIV)  ).GE.  0.0  )  GO  TO  530 
IFf  XU)  -  ABSI  HII. JPIV)  ) .  LT.  0.0  )  GO  TO  400 
480  HTEMP  =  HUPiV.I) 

Hi  IP  IV. I)  =  COS  I  NE  ’♦HTEMP  +  SI  NE*H(  JP  l  V,  I ) 

IF  I  XIIPIV  j  -  ABSI"  HUPIV.fj  ).GE.  0.0  )  GO  TO  500 

490  XIIPIV)  =  ABSI  HI  IPIV.I )  )  _  ....  __ 

IQ  I  IP  I V)  =  I 

500  HIJPIV.I)  =  —SI NE+HTEMP  ♦  COS INE +H IJPI V . I ) 

IF  I  XfjPIV)  -  A B SI  HUPiV.I)  t.uE.O.O)  GO"  TO  530 
510  XIJP1V)  =  ABSI  HIJPIV.I)  ) 

IOI  JP IV)  =  "I 

530  CONTINUE _ 

C  TEST  FOR  COMPUTATION  OF  EIGENVECTORS 

IFI IEGEN.NE.O)  GO  TO  40 
540  DO  550"  f  =  "i.N" 

HTEMP  =  UII. IPIV) 

II n.  IPi  Vi  =  COS  I  NE  +HTEMP  +  SI  NE  *Ul  I  .  JPI V  ) 

550  UU.JPIV)  f  -SINE*HTEMP+COSI  NE+UII  .JPIV) 

GO  TO  40 
1002  CONTINUE 
1000  RETURN 
END 
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SUBROUTINE  MESH 
IMPLICIT  R  EAL*8(  A-H,0— Z  I 

COMMON /R  1/0X0,0X149)  ,  D  V,BETA  ,RO,C,CS  ,XINF 
COMMON/I  L/N,M,NL 
COMMON /XX/XD*X( 50) 

IZ=0 

X{ IZ )=-X  INF 
XIN  )=  X  INF 
DO  1  1*1,5 
L  X(  I  )  =  -  XI  NF  + 1  *3*  DO 
DO  2  1=6, 44 

2  X<  I)=-5.D0  +  <  1-5)*. 25D0 
DO  3  1=45,49 

3  X<  I)=5.00+<  1-45  )*3.D0 
DO  4  I  =  I Z » N  L 

4  DX {  I  1=  XC  HII-XII  » 

RETURN 

END 


MULT 


SUBROUTINE  M UL T  ( A ,  B ) 

IMPLICIT  REAL* 8( A— H ,0-Zl 
COMMON/I  1/N,M,N1 

DIMENSION  A{49,  I)  ,8(49,1)  ,  TEMP  (49) 

DO  1  J  =  I  ,N  l 

DO  2  1=1, N1 

SUM=O.DO 

DO  3  K=l,Nl 

3  SUM=  SUM+A (1,K)*B(K,J) 

2  TEMPI  I  )=SUM 
DO  1  1=  1  ,N  l 
1  B(  1,  J  )  =  TEMP(  I) 

RETURN 

END 
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SUBROUTINE  MULT2(A,B) 

IMPLICIT  REAL*  8(  A— H  »  0-Z  > 

COMMON/I l/N,MfNl 

DIMENSION  A(49,1),B{49*II  ,  TEMP  (49) 

DO  1  I- L« N  1 

DO  2  J=1»N1 

SUM=O.DO 

DO  3  K= l»N I 

3  SUM=SUM+A( I »K1*B( K, J) 

2  TEMP  (  J  )=SUM 
DO  I  J  =  I f N  1 
I  A I  I f  J )  =  TEMP ( J ) 

RETURN 

END 


SETTUP 


SUBROUTINE  SETTUP 
IMPLICIT  REAL*8(A-H,0-Z» 

COMMON  /R  1/DX0»DX(49)  ,  DY,BfcTA  ,RO  ,C  ,CS  ,XI  NF,  PI2  *Y,DY2 
COMMON /i I/N«MfNlrMlyK|IDLfID2tI2fN2 
COMMON/BROWN /GAM 

ThIS  SUBROUTINE  IS  FOP.  INITIALIZING  PROGRAM  CONSTANTS 


Y=1.D0 
X  INF=  20.00 
N=  50 
M=  16 
K=  3 

8ETA=4  «D0 

R0=  1 *D0 

GAM= L.DO 

101=49 

102=50 

IZ=0 

N  I=N-  I 

N2=N-2 

0Y=2.D0*Y/M 

DY2=2.D0*DY 

PI=3.14I5  5265  3589  7932  DO 

P I  2=  2.D0*P I 

CALL  MESH 

CALL  FFIX 

RETURN 

END 
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TP  I  X 


SUBROUTINE  TFIX 
IMPLICIT  RfcAL*8( A-H, U-Z ) 

COMMON  T ( 49  I 
CGMMON/I  I / N » M *  N I 

CUMMON/R  I,/DX0,DX(4S)  ,JY,  BETA,  PO,C,CS  ,X  INF 
COMMON /XX/XD,XC  5  0) 

00  2  1*1, NI 
2  T (  I )  =TFUNC (  I  ,X(  I)) 

RETURN 

END 


TFUNC 

DOUBLE  PRECISION  FUNCTION  TFUNC(i,X) 

IMPLICIT  RE  AL*8(  A-H,0-Z  ) 

CGMMON/C  ONTRL  /  I  C 
T  FUNC  =  DM  A  X  1C  X,  O.DO) 

GO  TO  (  1,2,3,41,  IC 

1  RETURN 

2  I F ( I . EQ • 26) TFUNC =0.00 
IF  (  1  .FG.27)  TFJNC  =  .25U0 
I F ( I .EG.28JTFUNC*. 5D0 

F  F  TUP N 

3  iF(  1  .EG.  26.0R.  i  .EQ.27)  TFUNC  =  O.DO 
IP( I . EG . 28 ) TFUNC*. 25D  0 

I F ( I .E0.29)TFUNC=. 5DU 
IF ( i .EG. 30) TFUNC*. 85D0 
IF! I.E0.3l)T  FUNC* 1 • 2D  0 
IF (I .EG. 32) TFUNC* L. 5500 
RETURN 

4  IFC I .EG. 26 .OR. I .E 0.2 7. OR. I .EQ.28J TFUNC=O.DO 
I F { I .EQ.29)TFUNC=.2D0 

I F (  I . EQ . 30 ) TFUNC*. 45D  0 
IF  (  I  .EG.  3DTFUNC*.  7D0 
I F ( I.E0.32)TFUNC=. 97500 
IFC I . EG. 33) TFUNC* 1. 3D C 
IFC I .EG. 34) TFUNC* 1.60 C 
I  F  (  I.fq  .35  )TF(JNC  =  1.S5U0 
IFC I ,E0.36)TFUNC=2. 375D0 
IF ( 1 .EG. 37) TFUNC=2 • 85D0 
IFC  I  .EO  •  38  )  TFUNC*  3.  2  5D0 
RETURN 
END 
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NOMENCLATURE 

C  Airfoil  chord  length 

Cl  Lift  coefficient 

h  Semiheight  of  tunnel 

R{x)  Porosity  parameter 

5  Airfoil  surface  area 

U  Free- stream  velocity 

x,  y  Normalized  Cartesian  coordinates 

X,  Y  Cartesian  coordinates  (Fig.  1) 

T  |S/R(x) 

3  Compressibility  parameter 

7  r  Vortex  strength 

6  Lift  interference  factor,  Eq.  (13) 

6x,  6y  Finite  spacing  in  x  and  y  direction 

®  Perturbation  velocity  potential 

< j>  Interference  velocity  potential 

0  Model  velocity  potential 
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